source("00_util_scripts/mod_seurat.R")
source("00_util_scripts/mod_bplot.R")

proj.nm <- 'mission/SLE_TRPM2_MfMo/'

sprintf('%x',c(72,95,252))
kn.pal2 <- c('#d44','#485ffc')

# TRPM2 in all PBMC ----------
sobj <- read_zstd_rds('~/append-ssd/alaria2/GSE174188.perez.SLE.zst.rds')

sobj |> write_zstd_rds('~/append-ssd/alaria2/GSE174188.perez.SLE.zst.rds')

mex.perez <- sobj |> GetAssayData(layer = 'counts')

mex.perez |>
  DropletUtils::write10xCounts(path = '~/append-ssd/GSE174188.perez.SLE.h5',
                               x = _)

sobj <- sobj |>
  mutate(Status = fct_relevel(Status, 'Healthy', 'Treated', 'Managed'),
         manual.main = ifelse(str_detect(cov, 'dendritic'), 'Dendritic cells',
                              cov))

sobj |> DimPlot(group.by = 'manual.main', cols = 'Paired') +
  ggtitle('PBMC') +
  theme_jpub + NoLegend()

publish_pdf('mission/SLE_TRPM2_MfMo/perez.umap.pdf', height = 60, width = 57)

g0 <- sobj |> FeaturePlot('TRPM2', cols = c('lightgrey', 'red'), order = T,
                    split.by = 'SLE_status')

g0 + plot_annotation(title = 'TRPM2 expression') &
  theme_jpub(theme_classic) & NoLegend()

m2_bycell <- sobj |>
  get_abundance_sc_wide('TRPM2')

sobj$TRPM2 <- m2_bycell$TRPM2

sobj |>
  ggplot(aes(UMAP_1, UMAP_2, weight = TRPM2)) +
  geom_bin2d(bins = 256) +
  facet_wrap(~SLE_status) +
  scale_fill_gradient(low = 'lightgrey', high = 'red') +
  labs(title = 'TRPM2 expression') +
  theme_jpub(theme_classic) +
  NoLegend()

publish_pdf('sle.pbmc.trpm2.featureplot.pdf', width = 90)

m2_by_status <- sobj |>
  select(.cell, Status, TRPM2, UMAP_1, UMAP_2)
  
set.seed(42)

m2_by_status |>
  slice_sample(n = 10000, by = Status)

sobj |>
  mutate(Status = fct_relevel(Status, 'Healthy', 'Treated', 'Managed')) |>
  ggplot(aes(UMAP_1, UMAP_2, weight = TRPM2)) +
  geom_bin2d(bins = 256) +
  facet_wrap(~Status) +
  scale_fill_gradient(low = 'lightgrey', high = 'red') +
  labs(title = 'TRPM2 expression') +
  theme_jpub(theme_classic) +
  NoLegend()

publish_pdf('sle.pbmc.trpm2.featureplot.pdf', width = 100, height = 100)

g0 |> write_rds('mission/SLE_TRPM2_MfMo/results/pbmc.m2.featureplot.rds')

sobj |>
  DotPlot2d('TRPM2', Status, manual.main)

g1 <- last_plot() 

sle_pbmc_m2 |>
  filter(group.x %in% c('Healthy','Flare')) |>
  mutate(group.y = fct_reorder(group.y, avg.exp.scaled),
         group.x = ifelse(group.x == 'Flare', 'SLE', 'HC')) |>
  BubblePlot(d2 = T) +
  labs(y = 'Cell type', x = 'Group',
       title = 'TRPM2 expression in SLE PBMC') +
  theme_jpub() +
  RotatedAxis()

publish_source_plot('sle.pbmc.celltype.m2.dotplot', height = 55, width = 60)

sle_pbmc_m2 <-
  read_csv('mission/SLE_TRPM2_MfMo/results/sle.pbmc.celltype.m2.dotplot.csv')

sle_pbmc_m2

sobj <- read_rds("mission/SLE_TRPM2_MfMo/data/perez.mono.rds")

sobj |> DotPlot(c('TRPM2','CD14','FCGR3A'), cols = 'RdBu')

sobj |>
  DotPlot2d('TRPM2', Status, cov)

submono_m2 <- last_plot() |>
  pluck('data')

sle_pbmc_m2 |>
  bind_rows(submono_m2) |>
  filter(group.x == 'Monocytes')

## which cell type upregulate M2 higher? -------
Idents(sobj) <- 'cg_cov'

celltype.list <- Idents(sobj) |> unique()

deg.slevhc.pbmc <- celltype.list |>
  map(\(x)FindMarkers(sobj, group.by = 'Status', subset.ident = x,
                      ident.1 = c('Managed','Treated','Flare')) |>
        as_tibble(rownames = 'gene') |> mutate(cluster = x),
      .progress = T) |>
  list_rbind()

deg.slevhc.pbmc |>
  filter(p_val_adj < .05) |>
  write_csv('mission/perez_2022_SLE/perez.pbmc.slevhc.deg.csv')

deg.slevhc.pbmc <- read_csv('mission/perez_2022_SLE/perez.pbmc.slevhc.deg.csv')

deg.slevhc.pbmc |>
  filter(gene == 'TRPM2')

deg.flavhc.pbmc <- celltype.list |>
  map(\(x)FindMarkers(sobj, group.by = 'Status', subset.ident = x,
                      ident.1 = c('Flare','Treated'),
                      ident.2 = c('Healthy','Managed'), features = 'TRPM2') |>
        as_tibble(rownames = 'gene') |> mutate(cluster = x),
      .progress = T) |>
  list_rbind()

deg.flavhc.pbmc |>
  filter(gene == 'TRPM2')

deg.flavhc.pbmc |>
  mutate(cluster = case_match(cluster,
                              'Progen' ~ 'CD34+ progenitor',
                              'Prolif' ~ 'Proliferating lymphocytes',
                              'cM' ~ 'TRPM2-hi monocytes',
                              .default = cluster) |>
           fct_reorder(avg_log2FC)) |>
  ggplot(aes(x = avg_log2FC, y = cluster, fill = -log10(p_val_adj),
             label = str_c('p=',signif(p_val_adj, 2)))) +
  geom_col() +
  geom_text(nudge_x = 0.1, hjust = 'left') +
  scale_fill_distiller(palette = 'Reds', direction = 1) +
  theme_pubr() +
  labs(title = 'Fold change of TRPM2 in PBMC cell types: Flare vs Healthy & Managed')

## PBMC cell type change among SLE status ---------
sobj@meta.data |>
  as_tibble(rownames = '.cell') |>
  write_csv('DE_cells/results/perez_sle_sobj_meta.csv')

perez.meta <- read_delim('DE_cells/results/perez_sle_sobj_meta.csv.gz')

pbmc.conf <- perez.meta |>
  calc_frac_conf_on_grouped_count(ind_cov_batch_cov, cg_cov)

indiv.status <- perez.meta |>
  distinct(ind_cov_batch_cov, Status, Sex, pop_cov, Age)

pbmc.conf |>
  filter(cg_cov == 'ncM') |>
  left_join(indiv.status) |>
  mutate(Status = fct_relevel(Status, 'Healthy', 'Managed', 'Treated')) |>
  ggplot(aes(Status, fraction*100, fill = Status)) +
  geom_violin() +
  geom_jitter(width = .1, height = 0) +
  #stat_summary(geom = 'crossbar', fun = 'mean', width = .5) +
  theme_pubr() +
  scale_fill_brewer(palette = 'RdYlBu', direction = -1) +
  labs(title = 'Proportion of TRPM2-hi monocytes in PBMC',
       x = 'Donor group', y = '% in PBMC')

g1 <- last_plot()

g1$data |> write_csv('mission/SLE_TRPM2_MfMo/fraction.M2hi.PBMC.csv')

pbmc.conf |>
  filter(cg_cov == 'ncM') |>
  left_join(indiv.status) |>
  select(-c(cg_cov, conf.high, conf.low)) |>
  write_csv('mission/SLE_TRPM2_MfMo/results/fraction.M2lo.PBMC.csv')

# M2 in monocyte --------
sobj |> write_rds('mission/SLE_TRPM2_MfMo/data/perez.mono.rds')

sobj <- read_rds("mission/SLE_TRPM2_MfMo/data/perez.mono.rds")

sobj@meta.data |> head() |>
  as_tibble()

paired.pal <- RColorBrewer::brewer.pal(name = 'Paired', n = 12)

paired.pal <- c(paired.pal, '#AEC7E8', '#FFD7BE', '#47627A')

sobj |>
  DimPlot(cols = paired.pal, label = T, label.box = T, repel = T,
          label.size = 2) +
  theme_jpub

publish_pdf('perez.mono.umap.pdf', width = 60, height = 50)

## determine M2 subsets in mo --------
### average expr of M2 --------
sobj |> DotPlot('TRPM2', dot.scale = 5) +
  scale_color_gradient2(low = kn.pal2[2], high = kn.pal2[1]) +
  labs(x = 'Gene', y = 'Monocyte clusters')

publish_source_plot('sle.pbmc.mo.m2.dotplot',
                    width = 80, height = 75)

#### TRPM2 & non/classical mono marker ---------
sobj |> DotPlot(c('TRPM2','CD14','FCGR3A')) +
  scale_color_gradient2(low = kn.pal2[2], high = kn.pal2[1])

m2.ncmono <- last_plot() |>
  pluck('data')

m2.ncmono |>
  ggplot(aes(features.plot, id, fill = avg.exp.scaled)) +
  geom_tile() +
  scale_fill_gradient2(high = kn.pal2[1], low = kn.pal2[2]) +
  theme_jpub

m2.ncmono |>
  mutate(Gene = features.plot, `Monocyte clusters` = id,
         `Z-score` = avg.exp.scaled) |>
  tidyplot(Gene, `Monocyte clusters`, color = `Z-score`, width = 30) |>
  add_heatmap(rotate_labels = 0) |>
  adjust_colors(new_colors = c(kn.pal2[2],'white',kn.pal2[1])) |>
  save_plot('mission/SLE_TRPM2_MfMo/figs/perez.mono.m2.noncls.heatmap.pdf')

m2.clus.expr <- last_plot() |>
  pluck('data') |>
  as_tibble() |>
  mutate(id = fct_reorder(id, avg.exp.scaled))

### HMGB2 expr -------
sobj |> 
  DotPlot(c('TRPM2','HMGB2')) +
  scale_color_gradient2(low = '#00f', high = '#a00')

hmgb2.m2expr <- last_plot() |>
  pluck('data') |>
  as_tibble() 

hmgb2.m2expr |>
  pivot_wider(names_from = features.plot, values_from = avg.exp,
              id_cols = id) |>
  mutate(subtype = ifelse(TRPM2 > 0.05, 'Classical Mono', 'Non-classical Mono')) |>
  ggplot(aes(TRPM2, HMGB2)) +
  geom_smooth(method = 'lm', linetype = 'dashed', se = FALSE, color = 'grey') +
  geom_point(aes(color = subtype)) +
  geom_text_repel(aes(label = id), size = 2) +
  theme_bw() +
  stat_cor(size = 2) +
  labs(title = 'SLE PBMC monocyte clusters') +
  theme_jpub

hmgb2.m2expr |>
  ggplot(aes(TRPM2, HMGB2)) +
  geom_smooth(method = 'lm', linetype = 'dashed', color = 'grey') +
  geom_point(aes(color = subtype)) +
  geom_text_repel(aes(label = id), size = 2) +
  theme_classic(base_size = 6, base_family = 'ArialMT') +
  stat_cor(size = 2) +
  labs(title = 'SLE PBMC monocyte clusters') +
  theme_jpub

publish_source_plot('sle.pbmc.mo.trpm2.hmgb2.correlation', width = 70)

hmgb2.m2expr <-
  read_csv('mission/SLE_TRPM2_MfMo/results/sle.pbmc.mo.trpm2.hmgb2.correlation.csv')

sum22expr <- hmgb2.m2expr |>
  summarise(sum.exp = sum(avg.exp), .by = id)

#### HMGB2 across status -------
sobj |>
  DotPlot2d('HMGB2', Status, seurat_clusters)

### define double-high mono -------
sobj <- sobj |>
  mutate(m2b2.fine = ifelse(seurat_clusters %in% c(5:8,13,14), 'non-DH',
                            'DH'))

sobj |> DimPlot(group.by = 'm2b2.fine') +
  ggtitle('TRPM2-HMGB2 double high monocyte in SLE')

m2b2.expfc <- sobj |>
  FindMarkersAcrossVar(split.by = 'm2b2.fine', group.by = 'Status',
                       ident.1 = 'Healthy',
                       features = c('TRPM2','HMGB2'))

m2b2.expfc

sobj |>
  filter(m2b2.fine == 'DH') |>
  bill.violin('TRPM2', group.by = Status)

last_plot() |>
  pluck('data') |>
  select(3:4) |>
  write_source_csv('double.high.m2.expr.sle')

sobj |>
  filter(m2b2.fine == 'DH') |>
  bill.violin('HMGB2', group.by = Status)

last_plot() |>
  pluck('data') |>
  select(3:4) |>
  write_source_csv('double.high.hmgb2.expr.sle')

m2b2.conf <- sobj |>
  calc_frac_conf_on_grouped_count(ind_cov_batch_cov, m2b2.fine)

sample_status <- sobj |> distinct(ind_cov_batch_cov, Status)

m2b2.conf |>
  left_join(sample_status) |>
  ggplot(aes(Status, fraction)) +
  geom_boxplot() +
  facet_wrap(~m2b2.fine)

m2b2.conf |>
  left_join(sample_status) |>
  mutate(group = ifelse(Status == 'Healthy', 'HC', 'SLE')) |>
  filter(m2b2.fine == 'DH') |>
  tidyplot(group, fraction, color = group) |>
  add_data_points_beeswarm() |>
  add_mean_dash(color = 'black') |>
  add_test_pvalue(hide_info = T)

### CCR2 expr -------
sobj |> 
  DotPlot2d(c('TRPM2','CCR2'), group.x = seurat_clusters, group.y = SLE_status) +
  scale_color_gradient2(low = '#00f', high = '#a00')

ccr2.m2expr <- last_plot() |>
  pluck('data') |>
  as_tibble() 

ccr2.m2expr |>
  filter(group.y == 'SLE') |>
  pivot_wider(names_from = features.plot, values_from = avg.exp,
              id_cols = group.x) |>
  mutate(subtype = ifelse(TRPM2 > 0.04, 'Classical Mono', 'Non-classical Mono')) |>
  ggplot(aes(TRPM2, CCR2)) +
  geom_smooth(method = 'lm', linetype = 'dashed', se = FALSE, color = 'grey') +
  geom_point(aes(color = subtype)) +
  geom_text_repel(aes(label = group.x), size = 2) +
  theme_bw() +
  stat_cor(size = 2) +
  labs(title = 'SLE PBMC monocyte clusters') +
  theme_jpub

publish_source_plot('sle.pbmc.mo.trpm2.ccr2.correlation', width = 70)

ccr2.m2expr |>
  filter(group.y != 'SLE') |>
  pivot_wider(names_from = features.plot, values_from = avg.exp,
              id_cols = group.x) |>
  mutate(subtype = ifelse(TRPM2 > 0.04, 'Classical Mono', 'Non-classical Mono')) |>
  ggplot(aes(TRPM2, CCR2)) +
  geom_smooth(method = 'lm', linetype = 'dashed', se = FALSE, color = 'grey') +
  geom_point(aes(color = subtype)) +
  geom_text_repel(aes(label = group.x), size = 2) +
  theme_bw() +
  stat_cor(size = 2) +
  labs(title = 'HC PBMC monocyte clusters') +
  theme_jpub

publish_source_plot('HC.pbmc.mo.trpm2.ccr2.correlation', width = 70)

### CD14 expr -------
sobj |> 
  DotPlot2d(c('TRPM2','CD14'), group.x = seurat_clusters, group.y = SLE_status) +
  scale_color_gradient2(low = '#00f', high = '#a00')

cd14.m2expr <- last_plot() |>
  pluck('data') |>
  as_tibble() 

cd14.m2expr |>
  filter(group.y == 'SLE') |>
  pivot_wider(names_from = features.plot, values_from = avg.exp,
              id_cols = group.x) |>
  mutate(subtype = ifelse(TRPM2 > 0.04, 'Classical Mono', 'Non-classical Mono')) |>
  ggplot(aes(TRPM2, CD14)) +
  geom_smooth(method = 'lm', linetype = 'dashed', se = FALSE, color = 'grey') +
  geom_point(aes(color = subtype)) +
  geom_text_repel(aes(label = group.x), size = 2) +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  stat_cor(size = 2) +
  labs(title = 'SLE PBMC monocyte clusters') +
  theme_jpub

publish_source_plot('sle.pbmc.mo.trpm2.cd14.correlation', width = 70)

### M2 DEG SLE vs HC -------
m2.slevhc.leiden <-
sobj |> FindMarkersAcrossVar(split.by = 'seurat_clusters', group.by = 'Status',
                             ident.1 = 'Healthy', features = 'TRPM2',
                             logfc.threshold = 0, min.pct = 0)
-1:1 |>
  set_names() |>
  map(safely(\(x)tibble(x = log(x)))) |>
  list_transpose(simplify = T)

### cell fraction change --------
status.mo.conf <- sobj |>
  calc_frac_conf_on_grouped_count(Status, seurat_clusters)

pointsize <- AutoPointSize(sobj)

sobj |>
  left_join(status.mo.conf[c('Status','fraction','seurat_clusters')]) |>
  ggplot(aes(umap_1, umap_2, color = fraction)) +
  geom_point(size = pointsize) +
  facet_wrap(~Status)

status.mo.conf.umap <- last_plot() |>
  pluck('data')

status.mo.conf.umap |>
  ggplot(aes(umap_1, umap_2, color = fraction)) +
  geom_point(size = pointsize) +
  facet_wrap(~Status) +
  scale_color_distiller(palette = 'OrRd', direction = 1)

## M2 expr HC vs SLE --------
sobj |>
  bill.violin(features = 'TRPM2', group.by = Status)

sobj |>
  FindMarkersAcrossVar(split.by = 'cov', group.by = 'Status',
                       ident.1 = 'Healthy', features = 'TRPM2')

## M2 expr Healthy+Managed vs Flare+Treated --------
m2.sample <- sobj |>
  AverageExpression(features = 'TRPM2', assays = 'RNA',
                    group.by = c('ind_cov_batch_cov','Status'))

m2.cd14hl.sample <- sobj |>
  AverageExpression(features = 'TRPM2', assays = 'RNA',
                    group.by = c('ind_cov_batch_cov', 'cov'))

### pseudobulk ---------
mono.pbulk <- sobj |>
  AggregateExpression(return.seurat = T,
                      group.by = c('ind_cov_batch_cov','Status'))

res.pbulk <- mono.pbulk |>
  mutate(Status = str_extract(orig.ident, '[a-zA-Z]+$')) |>
  FindMarkers(group.by = 'Status', ident.1 = 'Healthy', test.use = 'DESeq2')

res.pbulk |>
  as_tibble(rownames = 'gene') |>
  filter(gene == 'TRPM2')

mono.pbulk |> get_abundance_sc_long('TRPM2')

m2.sample$RNA |>
  as_tibble() |>
  pivot_longer(everything()) |>
  mutate(Status = str_extract(name, '(Managed|Flare|Healthy|Treated)$')) |>
  write_csv('mission/SLE_TRPM2_MfMo/results/sle.pbmc.monocyte.trpm2.csv')

m2.sample <-
  read_csv('mission/SLE_TRPM2_MfMo/results/sle.pbmc.monocyte.trpm2.csv')

m2.sample |> arrange(name)

indiv.status |>
  mutate(ind_cov_batch_cov = str_replace_all(ind_cov_batch_cov, '_', '-')) |>
  separate(ind_cov_batch_cov, into = c('ind','batch'), sep = ':') |>
  mutate(name = str_c(batch, ind, Status, sep = '_')) |>
  arrange(name) |>
  inner_join(m2.sample) |>
  select(3:8) |>
  write_source_csv('sle.pbmc.mono.trpm2')

m2.cd14hl.sample <- m2.cd14hl.sample$RNA |>
  as_tibble() |>
  pivot_longer(everything()) |>
  mutate(cell_type = str_extract(name, '(nonc|class).+'),
         name = str_remove(name, '_(nonc|class).+')) |>
  pivot_wider(names_from = cell_type, values_from = value)

indiv.status |>
  mutate(ind_cov_batch_cov = str_replace_all(ind_cov_batch_cov, '_', '-')) |>
  separate(ind_cov_batch_cov, into = c('ind','batch'), sep = ':') |>
  mutate(name = str_c(batch, ind, sep = '_')) |>
  inner_join(m2.cd14hl.sample) |>
  select(3:9) |>
  write_source_csv('sle.pbmc.mono.subset.trpm2')

m2hi.hmvf <- m2.sample$RNA |>
  as_tibble() |>
  pivot_longer(everything()) |>
  mutate(celltype = ifelse(str_detect(name, 'nonclass'), 'TRPM2-low monocyte',
                           'TRPM2-high monocyte'),
         group = ifelse(str_detect(name, 'Healthy'), 'Healthy',
                        'SLE'))

m2hi.hmvf |>
  mutate(group = fct_relevel(group, 'Flare', after = 1)) |>
  ggplot(aes(group, value, color = group)) +
  geom_boxplot() +
  scale_color_hue(direction = -1) +
  theme_pubr() +
  facet_wrap(vars(celltype)) +
  stat_compare_means(comparisons = list(c('Healthy','SLE'))) +
  labs(y = 'Normalized expression', title = 'TRPM2 expression in PBMC monocytes')

m2hi.hmvf |>
  select(-name) |>
  pivot_wider(names_from = celltype, values_from = value,
              values_fn = list) |>
  unnest(-group) |>
  write_csv('mission/SLE_TRPM2_MfMo/results/sle.pbmc.monocyte.trpm2.boxplot.csv')

all.mono.marker <- sobj |> FindAllMarkers(only.pos = TRUE)

all.mono.marker |>
  filter(p_val_adj < 0.05) |>
  write_csv("mission/SLE_TRPM2_MfMo/all.mono.marker.csv")

sobj |> FindAllMarkers(features = "TRPM2")

publish_pdf("TRPM2.perez.mono.dotplot.pdf", 100, 100)

sobj <- sobj |> mutate(subgroup = str_c(Status, "_", seurat_clusters))

g1 <- sobj |> DotPlot("TRPM2", group.by = "subgroup")

g1$data |>
  separate(id, into = c("status", "cluster")) |>
  ggplot(aes(cluster, status, color = scale(avg.exp), size = pct.exp)) +
  geom_point() +
  scale_color_distiller(palette = "RdYlBu") +
  theme_pubr()

g1$data |>
  separate(id, into = c("status", "cluster")) |>
  filter(!(cluster %in% c("13", "14", "5", "7"))) |>
  ggplot(aes(cluster, status, color = scale(avg.exp), size = pct.exp)) +
  geom_point() +
  scale_color_distiller(palette = "RdYlBu") +
  theme_pubr()

m2.allmo.fc <- read_delim("mission/SLE_TRPM2_MfMo/M2.allmono.logfc.txt")
m2.allmo.fc |>
  ggplot(aes(as.character(cluster), avg_log2FC, fill = p_val_adj < .05)) +
  geom_col()

all.mono.marker |> filter(gene %in% c("IL6", "IFNG"))

sobj |> DimPlot(group.by = "cov")

sobj |> FindMarkers(features = "TRPM2",
                    group.by = "cov",
                    ident.1 = "classical monocytes")

sobj |> bill.violin("TRPM2", cov)

sobj |> VlnPlot("TRPM2", group.by = "cov")

status.all.conf <- sobj |>
  calc_frac_conf_on_grouped_count(group = Status, subtype = cov)

status.all.conf |>
  ggplot(aes(Status, fraction, fill = Status)) +
  geom_col() + facet_wrap(~cov, scales = "free_y") +
  theme_pubr()

sobj <- sobj |> mutate(ind_status = str_c(ind_cov, ":", Status))

ind_status_frac <- sobj |>
  calc_frac_conf_on_grouped_count(group = ind_status, subtype = cov) |>
  mutate(status = str_remove(ind_status, ".+:"))

ind_status_frac |> head()

ind_status_frac |>
  ggplot(aes(status, fraction, color = status)) +
  geom_boxplot() +
  geom_jitter(width = .2, height = 0) +
  facet_wrap(~cov, scales = "free_y") +
  theme_pubr() +
  stat_compare_means(comparison = list(c("Flare", "Healthy")))

sobj$SLE_status |> table()
sobj$Status |> table()

sobj |> FindMarkers(features = c("TRPM2", "IFNG"),
                    group.by = "SLE_status", ident.1 = "SLE")

sobj |> FindMarkers(features = c("TRPM2", "IFNG"),
                    group.by = "SLE_status", ident.1 = "SLE",
                    subset.ident = 1)

mono.leiden <- sobj |>
  FindMarkersAcrossVar(split.by = 'seurat_clusters', group.by = "Status",
                       ident.1 = "Flare", ident.2 = "Healthy")

flvhc.mono.fc |> filter(p_val_adj < .05) |>
  write_csv("mission/SLE_TRPM2_MfMo/perez.mono.flvhc.leiden.csv")

slevhc.mono.fc <-
  read_csv("mission/SLE_TRPM2_MfMo/perez.mono.slevhc.leiden.csv")

slevhc.mono.fc |>
  filter(gene %in% c("TRPM2", "IFNG", "IL6"))

sobj |> write_rds("mission/SLE_TRPM2_MfMo/perez.mono.rds")

cytok <- query_uniprot_keyword("KW-0202")

cytok

slevhc.mono.fc |>
  filter(gene %in% cytok$symbol)

flvhc.mono.fc |> filter(cluster == 1, gene %in% cytok$symbol) |>
  ggplot(aes(avg_log2FC, fct_reorder(gene, avg_log2FC))) +
  geom_col() + theme_pubr() +
  labs(title = "Cytokine in cluster 1 TRPM2-hi monocytes: SLE-flare vs HC", y = "Gene")

flvhc.mono.fc |> filter(cluster == 9, gene %in% cytok$symbol) |>
  ggplot(aes(avg_log2FC, fct_reorder(gene, avg_log2FC))) +
  geom_col() + theme_pubr() +
  labs(title = "Cytokine in cluster 9 TRPM2-hi monocytes: SLE-flare vs HC", y = "Gene")

slevhc.mono.fc |> filter(cluster == 11, gene %in% cytok$symbol) |>
  ggplot(aes(avg_log2FC, fct_reorder(gene, avg_log2FC))) +
  geom_col() + theme_pubr() +
  labs(title = "Cytokine in cluster 11 TRPM2-hi monocytes: SLE vs HC", y = "Gene")

sobj |> DotPlot("TRPM2", group.by = "Status", cols = "RdYlBu")

sobj |> bill.violin("TRPM2", Status)

sobj |> FindMarkers(features = "TRPM2", group.by = "Status",
                    ident.1 = "Flare", ident.2 = "Healthy")

sobj |> FindMarkers(features = "TRPM2", group.by = "Status",
                    ident.1 = "Flare", ident.2 = "Managed")

sobj |> FindMarkers(features = "TRPM2", group.by = "Status",
                    ident.1 = "Flare", ident.2 = "Treated")

## M2hi vs lo in SLE or in healthy ------------
sobj$cg_cov |> table()

stts.m2hvl.deg <- sobj |>
  FindMarkersAcrossVar(split.by = 'Status', group.by = 'cov',
                       ident.1 = 'classical monocytes')

stts.m2hvl.deg |>
  write_csv('mission/SLE_TRPM2_MfMo/results/SLE.status.m2hvl.deg.csv')

stts.m2hvl.deg <-
  read_csv('mission/SLE_TRPM2_MfMo/results/SLE.status.m2hvl.deg.csv')

flare.gsebp.tb <-
  read_csv('mission/SLE_TRPM2_MfMo/results/flare.m2hvl.gsebp.csv')

chemotax.gene <- flare.gsebp.tb |>
  filter(str_detect(Description, '^chemotaxis')) |>
  pull(core_enrichment) |>
  str_split_1('/')

inflam.gene <- flare.gsebp.tb |>
  filter(str_detect(Description, '^inflamma')) |>
  pull(core_enrichment) |>
  str_split_1('/')

### module score of inflammatory genes ----------
sobj <- sobj |>
  AddModuleScore(features = list(inflam.gene), name = 'inflam')

sobj |> VlnPlot('inflam1', pt.size = 0)

g1 <- last_plot()

inflam.module <- g1[[1]] |> pluck('data') |> as_tibble() |>
  summarise(inflam.score = mean(inflam1), .by = ident)

inflam.module |>
  write_csv('mission/SLE_TRPM2_MfMo/results/pbmc.mono.leiden.inflamma.module.csv')

inflam.module <-
  read_csv('mission/SLE_TRPM2_MfMo/results/pbmc.mono.leiden.inflamma.module.csv')

m2.clus.expr |>
  mutate(subtype = ifelse(avg.exp > .05, 'Classical Mono', 'Non-classical Mono')) |>
  left_join(inflam.module, join_by(id == ident)) |>
  ggplot(aes(avg.exp, inflam.score)) +
  geom_smooth(method = 'lm', linetype = 'dashed', se = FALSE, color = 'grey') +
  geom_point(aes(color = subtype)) +
  geom_text_repel(aes(label = id), size = 2) +
  theme_bw() +
  stat_cor(size = 2) +
  labs(x = 'Average expression of TRPM2', y = 'Inflammatory signature score',
       title = 'SLE PBMC monocyte clusters') +
  theme_jpub

sle.pbmc.mo.m2.inflam.score |>
  ggplot(aes(avg.exp, inflam.score)) +
  geom_smooth(method = 'lm', linetype = 'dashed', color = 'grey') +
  geom_point(aes(color = subtype)) +
  geom_text_repel(aes(label = id), size = 2) +
  theme_classic(base_size = 6, base_family = 'ArialMT') +
  stat_cor(size = 2) +
  labs(x = 'Average expression of TRPM2', y = 'Inflammatory signature score',
       title = 'SLE PBMC monocyte clusters') +
  theme_jpub

publish_source_plot('sle.pbmc.mo.m2.inflam.score.correlation', width = 70)

sle.pbmc.mo.m2.inflam.score <-
  read_csv('mission/SLE_TRPM2_MfMo/results/sle.pbmc.mo.m2.inflam.score.correlation.csv')

sum22expr |>
  mutate(subtype = ifelse(sum.exp > 2.5, 'Double-High Mono', 'Other Mono'),
         id = as.numeric(id)) |>
  left_join(inflam.module, join_by(id == ident)) |>
  ggplot(aes(sum.exp, inflam.score)) +
  geom_smooth(method = 'lm', linetype = 'dashed', se = FALSE, color = 'grey') +
  geom_point(aes(color = subtype)) +
  geom_text_repel(aes(label = id), size = 2) +
  theme_bw() +
  stat_cor(size = 2) +
  labs(x = 'Sum expression of TRPM2 + HMGB2', y = 'Inflammatory signature score',
       title = 'SLE PBMC monocyte clusters') +
  theme_jpub

publish_source_plot('sle.pbmc.mo.m2b2.inflam.score.correlation', width = 70)

### vlnplot of enriched genes ---------
stts.m2hvl.deg |>
  filter(gene %in% chemotax.gene, cluster == 'Flare') |>
  slice_max(avg_log2FC, n = 20, with_ties = F) |>
  mutate(p_val_adj = ifelse(p_val_adj == 0, 1e-300, p_val_adj),
         gene = fct_reorder(gene, avg_log2FC)) |>
  ggplot(aes(avg_log2FC, gene, fill = -log10(p_val_adj))) +
  geom_col() +
  scale_fill_distiller(palette = 'Reds', direction = 1) +
  theme_pubr(legend = 'right') +
  labs(title = 'Top enriched chemotaxis genes in flare SLE:\nTRPM2-hi vs TRPM2-lo') +
  theme_jpub

publish_source_plot('sle.pbmc.flare.chemotaxis.enrich.gene', width = 70)

stts.m2hvl.deg |>
  filter(gene %in% inflam.gene, cluster == 'Flare') |>
  slice_max(avg_log2FC, n = 20, with_ties = F) |>
  mutate(p_val_adj = ifelse(p_val_adj == 0, 1e-300, p_val_adj),
         gene = fct_reorder(gene, avg_log2FC)) |>
  ggplot(aes(avg_log2FC, gene, fill = -log10(p_val_adj))) +
  geom_col() +
  scale_fill_distiller(palette = 'Reds', direction = 1) +
  theme_pubr(legend = 'right') +
  labs(title = 'Top enriched inflammatory genes in flare SLE: TRPM2-hi vs TRPM2-lo')

top15.inflam <- stts.m2hvl.deg |>
  filter(gene %in% inflam.gene, cluster == 'Flare') |>
  slice_max(pct.1, n = 15, with_ties = F)

top15.inflam |>
  write_csv('mission/SLE_TRPM2_MfMo/results/top15.inflamm.gene.csv')

top15.inflam <-
  read.csv('mission/SLE_TRPM2_MfMo/results/top15.inflamm.gene.csv')

sobj <- sobj |>
  mutate(trpm2.type = ifelse(str_detect(cov, '^class'), 'M2-hi',
                             'M2-lo'))

sobj |>
  filter(Status == 'Flare') |>
  bill.violin(top15.inflam$gene, group.by = trpm2.type, facet.ncol = 5) +
  labs(x = 'Cell type', y = 'Normalized expression', fill = 'Cell type',
       title = 'Inflammatory response genes in SLE PBMC') +
  scale_fill_hue(labels = c('TRPM2-hi Mono', 'TRPM2-lo Mono')) +
  theme_jpub +
  theme(axis.text.x = element_blank())

publish_pdf('sle.pbmc.mo.inflamm.violin.pdf', width = 90)

## TRPM family ---------
trpmx <- sobj |> rownames() |> str_subset('TRPM\\d') |> str_sort()

trpmx.slevhc <- sobj |>
  FindMarkers(group.by = 'Status', ident.1 = 'Healthy',
              features = trpmx, min.pct = 0, logfc.threshold = 0)

trpmx.slevhc |>
  as_tibble(rownames = 'gene') |>
  mutate(logp = -log10(p_val_adj) * avg_log2FC / abs(avg_log2FC)) |>
  ggplot(aes(logp, gene, fill = avg_log2FC)) +
  geom_col() +
  geom_vline(xintercept = 0) +
  geom_vline(xintercept = -log10(.05), linetype = 'dashed') +
  scale_fill_gradient2(low = 'blue', high = 'red') +
  theme_pubr(legend = 'right') +
  labs(x = '-log10(adjusted P value)',
       title = 'TRPM family in monocytes from SLE PBMC',
       subtitle = 'GSE174188 (n=261)')

# heatmap ---------
m2hl.mrk.perez <- stts.m2hvl.deg |>
  mutate(direction = avg_log2FC > 0) |>
  filter(cluster == 'Flare', str_starts(gene, 'RP11-|MT-', negate = T)) |>
  slice_min(p_val_adj, n = 30, by = direction, with_ties = F) |>
  pull(gene)

ssobj <- sobj |>
  filter(Status == 'Flare')

g1 <- ssobj |>
  mutate(m2.type = ifelse(str_starts(cov, 'cla'),
                          'TRPM2-high', 'TRPM2-low')) |>
  ScaleData(features = m2hl.mrk.perez) |>
  DoHeatmap(m2hl.mrk.perez, group.by = 'm2.type', label = F, raster = F)

g1[[1]] + theme(axis.text.y = element_text(size = 2),
                text = element_text(size = 5))

publish_pdf('pbmc.m2.mono.heatmap.pdf', width = 90)

# surface marker of TRPM2-hi mono ------
stts.m2hvl.deg <-
  read_csv('mission/SLE_TRPM2_MfMo/results/SLE.status.m2hvl.deg.csv')

transmem <- query_uniprot_keyword('KW-0812')

gpi_anchor <- query_uniprot_keyword('KW-0336')

plsm_membr_go <- map_go_gene('GO:0005886')

trans_plsm_membr <- plsm_membr_go |>
  filter(SYMBOL %in% c(transmem$symbol, gpi_anchor$symbol))

trans_plsm_membr |>
  write_source_csv('plsm_membr_marker')

stts.m2hvl.mem <- stts.m2hvl.deg |>
  filter(avg_log2FC > 1, p_val_adj < .05) |>
  mutate(conserv_status = n(), .by = gene) |>
  filter(conserv_status == 4, gene %in% plsm_membr_go$SYMBOL) 

stts.m2hvl.mem |>
  write_csv('mission/SLE_TRPM2_MfMo/results/perez.mono.m2hvl.plasma.mem.csv')

stts.m2hvl.mem |>
  mutate(p_val_adj = ifelse(p_val_adj == 0, 1e-300, p_val_adj),
         gene = fct_reorder(gene, p_val_adj, .fun = mean, .desc = T)) |>
  slice_max(gene, n = 80) |>
  ggplot(aes(cluster, gene, fill = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point(shape = 21) +
  scale_fill_distiller(palette = 'Reds', direction = 1) +
  theme_bw() +
  labs(x = 'SLE status', subtitle = 'Sort by p-value',
       title = 'Top 20 surface protein of TRPM2-hi monocyte vs TRPM2-lo monocyte')

## correlation with TRPM2 --------
candid_mem <- stts.m2hvl.mem$gene |> unique()

sle_leiden_mean <- sobj |>
  filter(Status != 'Healthy') |>
  AverageExpression(features = plsm_membr_go$SYMBOL)

sle_leiden_mean <- sle_leiden_mean$RNA |>
  t() |>
  as_tibble(rownames = 'cluster')

sle_leiden_m2 <- sle_leiden_mean |>
  dplyr::select(cluster, TRPM2)

sle_m2_cor <- sle_leiden_mean |>
  pivot_longer(-1, names_to = 'gene') |>
  filter(gene %in% trans_plsm_membr$SYMBOL) |>
  left_join(sle_leiden_m2) |>
  summarize(correlation = cor(value, TRPM2),
            p.val = cor.test(value, TRPM2)$p.value,
            .by = gene) |>
  na.omit() |>
  mutate(rank = rank(correlation))
  
proj.nm <- 'mission/SLE_TRPM2_MfMo/'

sle_m2_cor |>
  arrange(correlation) |>
  write_source_csv('SLE.PBMC.membr.gene.cor.TRPM2.mono')

sle_m2_cor_min <- sle_m2_cor |>
  slice_min(correlation, n = 5)

sle_m2_cor_mm <- sle_m2_cor |>
  slice_max(correlation, n = 6) |>
  bind_rows(sle_m2_cor_min)

sle_m2_cor |>
  mutate(gene = fct_reorder(gene, correlation)) |>
  ggplot(aes(correlation, rank)) +
  geom_point(size = .01) +
  geom_point(data = sle_m2_cor_mm, size = 1, color = 'red') +
  geom_label_repel(data = sle_m2_cor_mm, aes(label = gene), size = 2,
                  nudge_y = ifelse(sle_m2_cor_mm$correlation > 0, -300, 300)) +
  theme_bw(base_size = 6) +
  labs(title = 'Membrane protein gene correlated with TRPM2 in SLE monocyte',
       fill = 'p value') +
  theme(plot.title.position = 'plot', legend.key.size = unit(4, 'mm'))

publish_source_plot('SLE.membrane.gene.cor.TRPM2.mono', width = 70)

sle_m2_cor |>
  mutate(gene = fct_reorder(gene, correlation)) |>
  slice_max(gene, n = 20) |>
  ggplot(aes(correlation, gene, fill = p.val)) +
  geom_point(shape = 21) +
  scale_fill_distiller(palette = 'Reds') +
  theme_bw(base_size = 6) +
  labs(title = 'Top 20 membrane protein gene correlated with TRPM2 in SLE monocyte',
       fill = 'p value') +
  theme(plot.title.position = 'plot', legend.key.size = unit(4, 'mm'))

publish_source_plot('top20.membrane.gene.cor.TRPM2.mono', width = 60)

sle_leiden_mean |>
  mutate(subtype = ifelse(TRPM2 > 0.04, 'Classical Mono', 'Non-classical Mono')) |>
  ggplot(aes(TRPM2, CD93)) +
  geom_smooth(method = 'lm', linetype = 'dashed', se = FALSE, color = 'grey') +
  geom_point(aes(color = subtype)) +
  geom_text_repel(aes(label = str_remove(cluster, 'g')), size = 2) +
  theme_bw() +
  stat_cor(size = 2) +
  labs(title = 'SLE PBMC monocyte clusters') +
  theme_jpub

publish_source_plot('sle.pbmc.mo.trpm2.cd93.correlation', width = 70)
